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Abstract 

This article concerns the application of bootstrap methodology to construct 
a likelihood-based confidence region for operating conditions associated with 
the maximum of a response surface constrained to a specified region. Unlike 
classical methods based on the stationary point, proper interpretation of this 
confidence region does not depend on unknown model parameters. In addition, 
the methodology does not require the assumption of normally distributed er- 
rors. The approach is demonstrated for concave-down and saddle system cases 
$_i ■ in two dimensions. Simulation studies were performed to assess the coverage 



probability of these regions. 

AMS 2000 subj classification: 62F25, 62F40, 62F30, 62J05. 

Key words: Stationary point; Kernel density estimator; Boundary kernel. 



*gibb. rd@pg.com, The Procter & Gamble Company, Mason, Ohio 45040-9452. 
t Correspondence: I-Li.Lu@boeing.com, Applied Statistics, Phantom Works, The Boeing Com- 
pany, P.O. Box 3707, MC 7L-22, Seattle, WA 98124-2207. 

■fwhcarter@vcu.edu, Department of Biostatistics, Virginia Commonwealth University, Richmond, 

VA 23298-0032. 

1 



2 Roger Gibb, I-Li Lu, and Walter Carter 

1 Introduction 



One of the reasons for using a response surface analysis is to determine the operating 
conditions which, without loss of generality, maximize a response y. Development 
of confidence regions for these operating conditions has been considered by several 
investigators. The purpose of this report is to demonstrate that likelihood-based 
bootstrap confidence region methodology is a powerful alternative which does not 
suffer from some of the limitations associated with existing approaches. 

Over the experimental region it is assumed that the relationship of y and the 
regressors x±, x 2 , . . . , x k can be expressed 



y = g(x,0) + e, 



(i) 



where g is an unknown continuous, differentiable function and e is a random source 
of variability not accounted for in g. Often g(x, 6) can be adequately approximated 
with a second-order polynomial, i.e. 



g{x) « ft + xf3 + x'Bx, 



(2) 



where (3 = (ft, ft, ■ ■ ■ , ft) and 



B = 



fti ft 2 /2 
ft 2 /2 ft 2 



Pik/2 
ft fe /2 



\ 



(3) 



\/V2 ft fc /2 ... (3 kk /2j 

Box and Hunter (1954) constructed a confidence region for the stationary point 
of a second-order response surface. The stationary point is given by 

B l (3 



(4) 
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the solution to the system of equations dy/dx = 0. The interpretation and relevance 
of the stationary point depends on the nature of the response surface which is deter- 
mined by B. Consider the following cases: 1) the eigenvalues of B are mixed in sign, 
2) the eigenvalues of B are all positive and 3) the eigenvalues of B are all negative. 

In cases 1 and 2 the stationary point is a saddle point and the location of minimum 
response, respectively, not the location of maximum response and, therefore, is not of 
interest for the purpose of this report. In case 3 the stationary point is the location 
of maximum response. In practice, the model parameters, including the elements 
of B, are unknown but can be estimated from the data. Since there is uncertainty 
associated with any estimator of B there is also uncertainty in assessing the nature 
of the true response surface and, therefore, the interpretation of Box and Hunter's 
confidence region. 

Another issue of practical importance is the data are observed over a treatment 
space of finite dimensions, i.e. the experimental region, and one is generally unwilling 
to extrapolate to areas outside this region. Peterson (1992) utilized a transforma- 
tion technique to develop a confidence region for the stationary point constrained 
to a specified region. However, as the approach is founded on the stationary point, 
the methodology suffers from the same interpretation difficulty as Box and Hunter's 
confidence region. 

In cases 1 and 2 the maximum response is undefined unless interest is restricted 
to a subset of the treatment space. Ridge analysis was developed by Hoerl (1959) 
and refined by Draper (1963) to estimate the stationary point subject to the con- 
straint x'x = r 2 . Under their approach the investigator must specify a value of the 
Lagrangian multiplier \i which, if chosen greater than the largest eigenvalue of B, 
facilitates determination of the location of the predicted constrained maximum. Sta- 
blein, Carter and Wampler (1983) constructed a confidence region for the constrained 
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stationary point conditional on the investigator's choice of ji. However, proper in- 
terpretation of their confidence region depends on whether the choice of \x is greater 
than the largest eigenvalue of B. Since the eigenvalues of B can only be estimated, 
there is uncertainty involved in the confidence region's interpretation. 

Let x cm be the operating conditions associated with maximum g(x, 6) subject to 
the constraint that x is within the experimental region. Denote x cm and as estima- 
tors of x cm and 6 respectively, then in practice, x cm can be calculated from g(x; 6) 
using a numerical optimization algorithm, such as the Nelder-Mead (1965) simplex. 
A favorable property of x cm is its interpretation does not depend on unknown model 
parameters, unlike the stationary point or constrained stationary point. A drawback, 
however, is that in some instances x cm is not a consistent estimator. Consider the 
case where g(x,0) is continuous and £j ~ N(0,o" 2 ). If the model is correct, 6 is a 
least-squares estimator and x cm is unique, then x cm can be shown to be consistent 
(Kendall and Stuart, 1979, Chapter 18). To demonstrate an instance where x cm is 
not unique and, therefore, x cm is not consistent, consider the k = 2 second-order 
response surface where fa — fa — and fa\ = P22 > 0. If the experimental region is 
the two dimensional direct product of the interval [—a, a], i.e. [—a, a] 2 , then the max- 
imum response within the experimental region occurs at (—a, —a), (a, —a), (—a, a) 
and (a, a). A formal assessment of whether x cm is unique could be made with tests 
of the appropriate hypotheses. For example, the likelihood of the symmetric model 
described earlier could be investigated by testing H : fa = fa = 0, fa\ = fai- If there 
is insufficient evidence to reject H further investigation is warranted. 

Construction of a confidence region for x cm using classical methods would require 
knowledge of the sampling distribution of x cm . As this distribution is unknown it 
is reasonable to explore the use of bootstrap confidence region methods that do not 
require its mathematical derivation. Unlike the confidence region methods described 
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earlier, the bootstrap approach does not require normality assumptions on the model 
errors. 



2 Likelihood-Based Bootstrap Confidence Regions 

Hall (1987, 1992) describes three methods for constructing likelihood-based bootstrap 
confidence regions for a A;-variate parameter vector £ given a sample of size n, namely, 
the percentile-t method, the ordinary percentile (hybrid) method, and the percentile 
method. In regards to which of these methods is appropriate when £ = x cm , there 
are several factors that merit consideration. First, the percentile-t method requires 
an accurate estimate for V, the asymptotic covariance matrix of n^(£ — £), where 
V is assumed to be positive definite. The analytic expression for V, an estimate of 
V, is unknown and accurate estimates using resampling techniques are unlikely to 
obtain in the small sample designed experiment setting. Second, and perhaps more 
importantly, only the percentile method preserves the range of £ in all cases. In 
particularly, for £ = x cm , only confidence regions under the percentile method are 
guaranteed to be bounded by the experimental region. Therefore, for the purposes 
of this report, attention is focused exclusively on the percentile method to construct 
likelihood-based bootstrap confidence regions for x cm . 

2.1 Percentile Method Algorithm 

An adaptation of Hall's (1987) likelihood-based bootstrap confidence region algorithm 
when £ = x cm is given below, followed by several notes regarding its implementation. 

I. For the response surface model y = X0 + e, determine using the method of 
least-squares. 
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2. Standardize the residual vector e with the elementwise operation 
e s = 6/ v /diag(/-X(X'X)- 1 X'). 

3. Compute b bootstrap estimates of x cm by following steps 3(a) - 3(d) a total of 
b times (see implementation notes for convenient choices of b). 

(a) Generate a bootstrap sample of e s , denoted by e*. 

(b) Calculate the corresponding bootstrap sample of the data with y* = X0 + 
e* 

(c) For the response surface model y* = X9* + e, determine using the 
method of least-squares. 

(d) Based on 6 , calculate x* cm . 

4. Fit a nonparametric density / to the x* cmi , i = 1, 2,. . . , b. 

5. The contour on / of smallest content that captures 100(1 — a)% of the b es- 
timates x* cm is a 100(1 — a)% likelihood-based confidence region for x cm using 
the percentile method. 

2.2 Implementation Notes 

Steps 1 — 3 are an application of bootstrapping residuals to generate bootstrap esti- 
mates of x cm . Bootstrapping residuals was proposed by Efron (1979) and is consid- 
ered appropriate when X is fixed and the model is correct with exchangeable errors. 
Standardization of the residual vector in step 2 ensures that the variance of each 
e s> i matches that of the unobservable E{. A simple bootstrap sample is generated 
by drawing a sample of size n with replacement from the elements of e s . The bal- 
anced bootstrap, proposed by Davison et al. (1986), was used in the examples and 
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simulation studies of section [2T41 to ensure that e = e. In step 4 the conditional dis- 
tribution of x cm is approximated nonparametrically and used as the basis for selecting 
the confidence region boundary. 

In practice, the level contour on / of smallest content that captures 100(1 — a)% 
of the b estimates x* cm was determined as follows: Choose b such that (1 — a)b is a 
non-negative integer and assume that f(x* cmi ), i = 1, 2,. . . , b, are unique. Let f a be 
the value of / corresponding to the level contour of smallest content on f(x) that 
captures exactly (1 — a)b of the x* cmi , i = 1, 2,. . . , b. To calculate f a consider that, 
by definition, f{x* cmi ) > f a for all x* cmi captured by the contour and f(x* cmi ) < f a 
for every other x* cmi . Therefore, f a is the (1 — a)b th element of the set with elements 
f(x* cmi ), i = 1, 2,. . . , b, sorted in descending order. In the event that all f(x* cmi ), 
i = 1, 2,. . . , b, are not unique, the contour f a (x) may capture more than (1 — a)b 
of the x* cmil i = 1, 2,. . . , b, in which case the confidence region is conservative. 
After identifying in this manner the confidence region boundary was graphically 
displayed using the GCONTOUR procedure of SAS®. 



2.3 Application to Second-Order Response Surfaces 

Likelihood-based confidence regions were constructed for x cm , where g(x,0) was a 
second-order polynomial (k = 2) for the following cases: 1) concave down with x cm 
located inside the experimental region and 2) saddle system with x cm located, neces- 
sarily, on a boundary of the experimental region. The experimental region was defined 
as the two dimensional direct product of the interval [—1.4, 1.4]. Model parameters 
and coordinates of x cm for both models are summarized in table [TJ 

As indicated in section I2.1[ a nonparametric estimate of the sampling distribution 
of x cm conditional on b bootstrap estimates of x cm serves as the basis for constructing 
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the confidence region boundary. A product kernel density estimator for this purpose 
is given by 

/(*H^£{n<^H}. o» 

where K is a univariate kernel and hi is the bandwidth for the zth coordinate (see, 
e.g., Silverman, 1986, Scott, 1992). Gasser and Miiller (1979) showed that if K is 
symmetric the bias of f(x) can be inflated near the boundaries. Miiller (1988) showed 
that uniform bias over the entire support of f{x) can be achieved if an asymmetric 
boundary corrected kernel is used in (jSJ). Gasser and Miiller (1979) and Jones (1993), 
among others, describe univariate boundary kernels designed for estimating densities 
with a single boundary. Densities with support over an interval may be estimated 
using two such kernels, one at each boundary, but this is appropriate only if the 
bandwidth is small relative to the interval length. Hart and Wehrly's (1992) linear 
boundary kernel was implemented in (jHJ), as this kernel is specifically designed to 
estimate densities with support on a interval, where the bandwidth need not be small 
relative the interval length. The biweight density was used as a basis for constructing 
this boundary kernel. 

Two methods of multivariate bandwidth selection were explored: Normal rule-of- 
thumb (see, e.g., Scott, 1992) and the plug-in approach of Wand and Jones (1994, 
pp 107-113). The sample standard deviation in each coordinate of x* cm was used to 
estimate the scale of the conditional distribution of x cm , as some measure of scale 
is required by both bandwidth selectors. The possibility exists for all bootstrap 
estimates for x cm to be located on the experimental region boundary, a situation that 
occured almost exclusively with the saddle system. To ensure nonzero scale estimates, 
the standard deviation for the ith coordinate was calculated from a modified dataset 
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given by 




otherwise, 



1.4 



(6) 



where 5 is a random [7(0,0.05) deviate. The effect of using the modified dataset 
to estimate standard deviation was negligible in the concave-down system, as most 
bootstrap estimates did not fall on a boundary. For the saddle system, using the 
modified dataset prevented nonzero bandwidth estimates and, thus, permitted use of 
the product kernel density estimator. 

2.4 Simulation Study Results 

It is of interest to compare the bootstrap confidence region performance under both 
methods of bandwidth selection, i.e. Normal rule-of-thumb and the plug-in approach. 
Data for 500 experiments were simulated by adding a N(0, 3 2 ) deviate to the 'true' 
response at operating conditions associated with a k = 2 rotatable central-composite 
design (5 center runs, n = 13). Confidence regions using both bandwidth selectors 
were constructed using 2000 bootstrap samples in each simulated experiment. Band- 
widths under the plug-in method were slightly larger, on average, than those of the 
Normal rule-of-thumb (see table [2]). This result is not unexpected, as simulation 
studies of Wand and Jones (1994) indicate a tendency for the multivariate plug-in 
method to oversmooth in some cases. 

A contour of the true concave-down response surface and a representative con- 
fidence region for x cm with b = 2000 are provided in figures [1] and El respectively. 
Bootstrap estimates for x cm are indicated with points in figure El with a solid line 
identifying the confidence region boundary. Note that the confidence region bound- 
ary is within the experimental region, as required. Had the hybrid or percentile-t 
methods been implemented this would not have been true in all cases. 
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The simulated experiments were arranged into five groups of 100 and the coverage 
probability calculated in each group for confidence coefficients ranging from 90—100%. 
The mean coverage probability, with error bars placed at one standard error from the 
mean in either direction, is plotted in figure [31 Note that there is little difference 
in coverage probability under the two approaches, a result expected from the close 
similarity in bandwidths seen in table [2J One approach for potentially improving the 
coverage probability is explored in the discussion section. 

A contour plot of the true saddle system and a representative confidence region 
for x cm are provided in figures H] and El respectively. Note that in this case x cm is 
located on the boundary of the experimental region. Since all bootstrap estimates for 
x cm are also on the boundary, visual identification of the confidence region is difficult. 
In the figure, two small triangles identify the location of lower and upper boundaries 
in the x<i coordinate. In the X\ coordinate the confidence region extends slightly away 
from the experimental region boundary. 

Figure [6] summarizes the coverage probability under both bandwidth selectors 
for the saddle system. As with the concave-down response surface, the coverage 
probability is statistically indistinguishable under the two approaches. The coverage 
probability is closer to the nominal level than for the concave-down response surface. 

Calculation of the Normal rule-of-thumb bandwidths requires less time and com- 
putational resources than bandwidths under Wand and Jones' plug-in method. Since 
the coverage probability was essentially the same under both methods, the Normal 
rule-of-thumb method was implemented in the simulation studies that follow. 

A simulation study was conducted to assess whether increasing the number of 
bootstrap samples would improve the coverage probability. Confidence regions were 
constructed from the same simulated data described earlier using b = 4000 and b = 
6000. In neither the saddle system nor concave-down case was the coverage probability 
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significantly effected by increasing b. 

To explore the effect of sample size on coverage probability, confidence regions for 
x cm were constructed under both second-order response surfaces for n = 13, n = 26 
and n = 208, corresponding to 1, 2 and 16 experimental replications of a rotatable 
central-composite design with 5 center runs. The results for the concave-down and 
saddle system with b = 2000 are summarized in figures [7] and [BJ respectively. In both 
cases the coverage probability approaches the nominal level with increasing sample 
size. 

3 Discussion 

Existing confidence regions for operating conditions associated with the maximum 
of a response surface, either unconstrained or constrained to a specified region, are 
based on the stationary point. These approaches are only applicable to second-order 
models under the assumption of normally distributed errors. Also, the interpretation 
of these confidence regions can be ambiguous, since this requires assessment of the 
unknown elements of the B matrix. 

In contrast, the interpretation of a likelihood-based bootstrap confidence region 
for x cm does not depend on the nature of the response surface, assuming x cm is 
unique. For example, if g(x,0) is a second-order polynomial the confidence region 
interpretation is the same whether g(x, 6) is concave-down, concave-up or a saddle 
system. In addition, the approach is not restricted to second-order models nor does 
it require assumptions on the model error distribution, except for exchangeability of 
the errors. In principle, the methodology is also applicable to models where g(x, 6) is 
nonlinear in 6, though Hjorth (1994, pp 190) indicates that in non-linear regression 
applications a direct analogue of standardized residuals is generally not available. 
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Simulation results from section [2^1 (see figure [7]) indicate the bootstrap confidence 
region coverage probability was less than nominal under the concave-down system 
(n = 13). Coverage probability was higher for the saddle system (n = 13), but still 
below the nominal level (see figure [8]). A possible reason that higher coverage proba- 
bility was observed for the saddle system is the estimated conditional distribution of 
x cm is essentially restricted to one dimension in this case, thereby reducing the 'curse 
of dimensionality'. Evidence that the coverage probability in both systems converges 
to the nominal level with increasing sample size, an indication of the asymptotic 
accuracy of the methodology, is apparent in figures [7] and El 

Loh (1987, 1991) proposed a bootstrap calibration technique to improve the ac- 
curacy of confidence sets. In brief, the approach consists of estimating the confidence 
coefficient associated with the desired coverage probability. For example, to achieve a 
coverage probability of 90% in the concave-down response surface case (n = 13), it is 
apparent in figure [7] that a confidence coefficient of approximately 99.5% is required. 
It is also apparent that in this case realization of coverage probabilities greater than 
~ 91% are not feasible with this approach, unless the sample size is increased. For 
the saddle system (n = 13) a confidence coefficient of approximately 95% is needed 
to achieve a coverage probability of 90% (see figure E]). 

The accuracy of bootstrap confidence regions for x cm using the percentile method 
is largely determined by how accurately the conditional distribution of x cm , f, is 
estimated. Under the kernel method implemented in section 2, choice of bandwidth 
has a critical impact on this accuracy. Hall (1987) indicated that bandwidths which 
are optimal in a global sense, such as minimization of the mean integrated squared 
error of /, can produce confidence regions that are "unduly bumpy". He noted that 
this can be attributed to the fact that the ratio of the variance to the squared bias 
of / is greater in the tails than in the central region of the distribution. Therefore, 
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the Normal rule-of-thumb and plug-in bandwidth selectors were investigated for their 
tendency to avoid undersmoothing. The k = 2 concave-down and saddle system 
cases explored in section 12.41 are examples where these methods provide reasonable 
bandwidth estimates. Jones (1993) followed a similar approach when he used a plug-in 
bandwidth selector in conjunction with a univariate boundary kernel. 

There are situations, however, in which bandwidth selection is more challenging. 
For example, if / is bimodal, use of the standard deviation to estimate the scale of 
/ can result in considerable oversmoothing. Janssen et al. (1995) developed more 
robust scale estimators for such cases. Consider, also, a situation where all bootstrap 
estimates for x cm are equally distributed on the two boundaries A = {x: x\ = 
lA,x 2 e (-1.4,1.4)} and B = {x: x\ G (—1.4, 1.4), x 2 = 1.4}. On boundary A 
bandwidth h\ should be much smaller than h 2 . However, on boundary B the opposite 
is true. In this variable kernel density estimator would seem more appropriate 

than the fixed bandwidth kernel estimator implemented in section \2A\ as this would 
allow for different levels of smoothing depending on the location of bootstrap estimates 
for x cm . 

Application of the bootstrap confidence region methodology was restricted to the 
k = 2 case where the experimental region was rectangular in shape. The bootstrap 
and kernel density estimation techniques described in section 2 are not limited to two 
dimensions. Wand and Jones' plug-in bandwidth selector can also be extended to 
higher dimensions, though with greater computational expense. In light of the close 
comparison of coverage probability under the Normal rule-of-thumb and plug-in meth- 
ods observed in section 12.41 rule-of-thumb methods may be adequate in many higher- 
dimensional applications. Graphical presentation of confidence regions for k > 3 is a 
challenging issue, though this is not unique to our application. Staniswalis, Messer 
and Finston (1993) describe a multivariate boundary kernel designed for regions of 
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arbitrary shape which obtains the correct order of bias over the entire region (see, 
also, Scott, 1992, pp 155). 

Throughout this article discussion has focused on maximizing a single response. 
However, researchers in many areas of application are faced with the problem of si- 
multaneously improving multiple responses that depend on a common set of control- 
lable variables. Since a single operating condition is rarely optimal for all responses, 
compromise must be incorporated into the estimation procedure. Desirability opti- 
mization methodology (Harrington, 1965, Derringer and Suich, 1980) addresses this 
problem and has been proven effective in a wide range of applications involving contin- 
uous responses. In brief, the approach consists of estimating the operating conditions 
that maximize D, a simultaneous measure of the desirability of all the responses. An 
arguable shortcoming of the methodology is it does not provide for estimation of the 
variability of these operating conditions. Bootstrap confidence region methodology is 
defined in sufficiently general terms to allow for construction of a confidence region 
for the operating condition that maximizes D constrained to the experimental region, 
under the assumption that this parameter is unique. 

Of practical interest is the computational resources and time required to construct 
a likelihood-based bootstrap confidence region. The programs that generated the 
results of section 12.41 were written in SAS® and run on a 300 megahertz Pentium® II 
desktop computer with 128 megabytes of RAM. Approximately 2 minutes is required 
to generate a single confidence region with b = 2000, of which about 25 seconds are 
used for steps 1 — 3 with the balance of the time taken in steps 4 and 5. Wand (1994) 
describes methods of accelerating computations required for multivariate bandwidth 
optimization and kernel density estimation. These methods were not implemented 
due to technical limitations of SAS®. It is expected that steps 4 and 5 would require 
considerably less time under their optimized procedures. 
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Table 1: True second-order model parameters and associated x. 



Response 
















Surface 


Po 


A 


ft 


012 


0n 


022 




concave down 


86.850 


5.242 


4.778 


-0.775 


-2.781 


-2.524 


(0.828,0.819) 


saddle 


90.259 


-6.425 


1.244 


-0.775 


2.781 


-2.524 


(-1.4,0.462) 



Table 2: Mean and standard error (in parenthesis) of bandwidths estimated from 500 
simulated response surfaces using the Normal rule-of-thumb and plug-in methods. 





Concave-down 


Saddle system 


Bandwidth selector 


h 


h 2 


hi 


h 2 


Normal rule-of-thumb 


0.196 


0.213 


0.011 


0.261 




(0.0033) 


(0.0037) 


(6.9 x 10~ 6 ) 


(0.0049) 


Plug-in method 


0.214 


0.233 


0.013 


0.313 




(0.0034) 


(0.0038) 


(13 x 10- 6 ) 


(0.0060) 
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Figure 2: A 90% confidence region for x cm where the true response surface is concave- 
down and x identifies x cm . 
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Figure 3: Comparison of coverage probability under the Normal rule-of-thumb and 
plug-in bandwidth selectors. The true response surface is concave-down. 
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Figure 4: Contour plot of the true saddle system response surface, where x identifies 
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Figure 5: A 90% confidence region for x cm where the true response surface is saddle 
system and x identifies x cm . 
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Figure 6: Comparison of coverage probability under the Normal rule-of-thumb and 
plug-in bandwidth selectors. The true response surface is a saddle system. 
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Figure 7: Comparison of coverage probability for three different sample sizes where 
the response surface is concave-down and b = 2000 with Normal-rule-of-thumb band- 
widths. 
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Figure 8: Comparison of coverage probability for three different sample sizes where 
the response surface is a saddle system and b = 2000 with Normal-rule-of-thumb 
bandwidths. 



